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BACKGROUND 

Field of the Invention 

The present invention relates to performing arithmetic operations on 
interval operands within a computer system. More specifically, the present 
invention relates to a method and an apparatus that uses interval arithmetic to 
bound the solution set of a system of linear eqxxations. 

Related Art 

Rapid advances in computing technology make it possible to perform 
trillions of computational operations each second. This tremendous 
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computational speed makes it practical to perform computationally intensive tasks 
as diverse as predicting the weather and optimizing the design of an aircraft 
engine. Such computational tasks are typically performed using machine- 
representable floating-point numbers to approximate values of real numbers. (For 
5 example, see the Institute of Electrical and Electronics Engineers (IEEE) standard 
754 for binary floating-point numbers.) 

In spite of their limitations, floating-point numbers are generally used to 
perform most computational tasks. 

One limitation is that machine-representable floating-point numbers have a 

10 fixed-size word length, which limits their accuracy. Note that a floating-point 
number is typically encoded using a 32, 64 or 128-bit binary number, which 
means that there are only 2^^, 2^^ or 2^^^ possible symbols that can be used to 
specify a floating-point number. Hence, most real number values can only be 
approximated with a corresponding floating-point number. This creates 

1 5 estimation errors that can be magnified through even a few computations, thereby 
adversely affecting the accuracy of a computation. 

A related limitation is that floating-point numbers contain no information 
about their accuracy. Most measured data values include some amount of error 
that arises from the measurement process itself This error can often be quantified 

20 as an accuracy parameter, which can subsequently be used to determine the 

accuracy of a computation. However, floating-point numbers are not designed to 
keep track of accuracy mformation, whether fi^om input data measurement errors 
or machine rounding errors. Hence, it is not possible to determine the accuracy of 
a computation by merely examining the floating-point number that results from 

25 the computation. 

Interval arithmetic has been developed to solve the above-described 
problems. Interval arithmetic represents numbers as intervals specified by a first 
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(left) endpoint and a second (right) endpoint. For example, the interval [a, b], 
where a < 6, is a closed, bounded subset of the real numbers, R, which includes a 
and b as well as all real numbers between amdb. Arithmetic operations on 
interval operands (interval arithmetic) are defined so that interval results always 
5 contain the entire set of possible values. The result is a mathematical system for 
rigorously bounding numerical errors from all sources, including measurement 
data errors, machine rounding errors and their interactions. (Note that the first 
endpoint normally contains the "infimum", which is the largest number that is less 
H than or equal to each of a given set of real numbers. Similarly, the second 

p 1 0 endpoint normally contains the "supremum", which is the smallest number that is 

fl greater than or equal to each of the given set of real numbers.) 

=^ One commonly performed computational operation is to find the sohition 

i of a system of interval linear equations. What is needed is a method and an 

apparatus that uses interval arithmetic to efficiently compute narrow bounds on 
1 5 the solution set of a system of linear equations. 

Ui 

I SUMMARY 

One embodiment of the present invention provides a system that bounds 
the solution set of a system of linear equations Ax = b, wherein A is an interval 

20 matrix and b is an interval vector. During operation, the system preconditions the 
set of linear equations Ax = b by multiplying through by a matrix B to produce a 
preconditioned set of linear equations Mqx = r, wherein Mo = BA and r = Bb. 
Next, the system widens the matrix Mo to produce a widened matrix M, wherein 
the midpoints of the elements of M form the identity matrix. Finally, the system 

25 uses M and r to compute the hull h of the system Mx = r, which bounds the 
solution set of the system Mqx = r. 
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In a variation on this embodiment, the system computes the matrix B by 
computing an approximate center Ac of the matrix A, and then forming B by 
computing the approximate inverse of Ac, B = (Ac)"^ 

In a variation on this embodiment, the system additionally assures that 
5 supin) > 0 by changing the sign of and if necessary. 

In a variation on this embodiment, the system uses M and r to compute the 
hull h by forming P as an inverse of the left endpoint of M. The system also 
forms Cj = l/(2Pa - 1) for / = 1 , . . . , w and forms z^={inj{r)+sup{r))Pu - e^Ysup{r\ 
wherein e^^is a unit vector in which the z-th element is 1 and other elements are 0. 
10 The system then forms h by: setting infih^ = c,z, if > 0; setting infih) = z, \izi 
< 0; and setting sup(h) = Vsup(r), 

In a variation on this embodiment, the system determines whether or not 
M is regular. If the inverse of inf(M) exists and is denoted by P, and if 
/w/(M„) > 0 for all z, then M, Mo and A are all regular if and only if P > I. If not, 
1 5 the system terminates the process of computing the hull h. 

BRIEF DESCRIPTION OF THE FIGURES 

FIG. 1 illustrates a computer system in accordance with an embodiment of 
the present invention. 
20 FIG. 2 illustrates the process of compiling and using code for interval 

computations in accordance with an embodiment of the present invention. 

FIG. 3 illustrates an arithmetic unit for interval computations in 
accordance with an embodiment of the present invention. 

FIG. 4 is a flow chart illustrating the process of performing an interval 
25 computation in accordance with an embodiment of the present invention. 

FIG. 5 illustrates four different interval operations in accordance with an 
embodiment of the present invention. 
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FIG. 6 illustrates the process of bounding the solution set of a system of 
linear equations in accordance with an embodiment of the present invention. 



DETAILED DESCRIPTION 

5 The following description is presented to enable any person skilled in the 

art to make and use the invention, and is provided in the context of a particular 
application and its requirements. Various modifications to the disclosed 
embodiments will be readily apparent to those skilled in the art, and the general 
principles defined herein may be applied to other embodiments and applications 

1 0 without departing from the spirit and scope of the present invention. Thus, the 
present invention is not limited to the embodiments shown, but is to be accorded 
the widest scope consistent with the principles and features disclosed herein. 

The data structures and code described in this detailed description are 
typically stored on a computer readable storage medium, which may be any device 

15 or medium that can store code and/or data for use by a computer system. This 
includes, but is not limited to, magnetic and optical storage devices such as disk 
drives, magnetic tape, CDs (compact discs) and DVDs (digital versatile discs or 
digital video discs), and computer instruction signals embodied in a transmission 
medium (with or without a carrier wave upon which the signals are modulated), 

20 For example, the transmission medium may include a communications network, 
such as the Intemet. 



Computer System 

FIG. 1 illustrates a computer system 100 in accordance with an 
25 embodiment of the present invention. As illustrated in FIG. 1, computer system 
100 includes processor 102, which is coupled to a memory 112 and a to peripheral 
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bus 1 10 through bridge 106. Bridge 106 can generally include any type of 
circuitry for coupling components of computer system 100 together. 

Processor 102 can include any type of processor, including, but not limited 
to, a microprocessor, a mainframe computer, a digital signal processor, a personal 
organizer, a device controller and a computational engine within an appliance. 
Processor 102 includes an arithmetic unit 104, which is capable of performing 
computational operations using floating-point numbers. 

Processor 102 communicates with storage device 108 through bridge 106 
and peripheral bus 110. Storage device 108 can include any type of non-volatile 
storage device that can be coupled to a computer system. This includes, but is not 
limited to, magnetic, optical, and magneto-optical storage devices, as well as 
storage devices based on flash memory and/or battery-backed up memory. 

Processor 102 communicates with memory 1 12 through bridge 106. 
Memory 1 12 can include any type of memory that can store code and data for 
execution by processor 102. As illustrated in FIG. 1, memory 1 12 contains 
computational code for intervals 1 14. Computational code 114 contains 
instructions for the interval operations to be performed on individual operands, or 
interval values 115, which are also stored within memory 112. This 
computational code 114 and these interval values 1 15 are described in more detail 
below with reference to FIGs. 2-5, 

Note that although the present invention is described in the context of 
computer system 100 illustrated in FIG. 1, the present invention can generally 
operate on any type of computing device that can perform computations involving 
floating-point numbers. Hence, the present invention is not limited to the 
computer system 100 illustrated in FIG. 1. 
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Compiling and Using Interval Code 

FIG. 2 illustrates the process of compiling and using code for interval 
computations in accordance with an embodiment of the present invention. The 
system starts with source code 202, which specifies a number of computational 
5 operations involving intervals. Source code 202 passes through compiler 204, 
which converts source code 202 into executable code form 206 for interval 
computations. Processor 102 retrieves executable code 206 and uses it to control 
the operation of arithmetic unit 104. 

Processor 102 also retrieves interval values 115 from memory 112 and 
10 passes these interval values 115 through arithmetic unit 104 to produce results 
212. Results 212 can also include interval values. 

Note that the term "compilation" as used in this specification is to be 
construed broadly to include pre-compilation and just-in-time compilation, as well 
as use of an interpreter that interprets instructions at run-time. Hence, the term 
15 "compiler" as used in the specification and the claims refers to pre-compilers, 
just-in-time compilers and interpreters. 

Arithmetic Unit for Interv als 

FIG. 3 illustrates arithmetic unit 104 for interval computations in more 
20 detail accordance with an embodiment of the present invention. Details regarding 
the construction of such an arithmetic unit are well knovra in the art. For 
example, see U.S. Patent Nos. 5,687,106 and 6,044,454, which are hereby 
incorporated by reference in order to provide details on the construction of such 
an arithmetic unit. Arithmetic unit 104 receives intervals 302 and 312 as inputs 
25 and produces interval 322 as an output. 

In the embodiment illustrated in FIG. 3, interval 302 includes a first 
floating-point number 304 representing a first endpoint of mterval 302, and a 
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second floating-point number 306 representing a second endpoint of interval 302. 
Similarly, interval 312 includes a first floating-point number 314 representing a 
first endpoint of interval 312, and a second floating-point number 316 
representing a second endpoint of interval 312. Also, the resulting interval 322 
5 includes a first floating-point number 324 representing a first endpoint of interval 
322, and a second floating-point number 326 representing a second endpoint of 
interval 322. 

Note that arithmetic unit 104 includes circuitry for performing the interval 
operations that are outlined in FIG, 5. This circuitry enables the interval 
10 operations to be performed efficiently. 

How^ever, note that the present invention can also be applied to computing 
devices that do not include special-purpose hardware for performing interval 
operations. In such computing devices, compiler 204 converts interval operations 
into a executable code that can be executed using standard computational 
1 5 hardware that is not specially designed for interval operations. 

If FIG. 4 is a flow chart illustrating the process of performing an interval 

tn 

=£i computation in accordance with an embodiment of the present invention. The 

system starts by receiving a representation of an interval, such as first floating- 
point number 304 and second floating-point number 306 (step 402). Next, the 
20 system performs an arithmetic operation using the representation of the interval to 
produce a result (step 404). The possibilities for this arithmetic operation are 
described in more detail below with reference to FIG. 5. 



Interval Operations 

25 FIG. 5 illustrates four different interval operations in accordance with an 

embodiment of the present invention. These interval operations operate on the 
intervals X and 7. The interval X includes two endpoints, 

8 

Attorney Docket No. SIJN-P643 1 -SPL Inventors: Walster et al. 

ARPH \SUN MICROSYSTEMS\SUN-P6431-SPL\SUN-P6431-SPL APPLICATION DOC 



X denotes the lower bound of X, and 
X denotes the upper bound of X. 

The interval Xis a closed subset of the extended (including and +^) real 
numbers jR* (see line 1 of FIG. 5). Similarly the interval Falso has two endpoints 
and is a closed subset of the extended real numbers 7?* (see line 2 of FIG. 5). 

Note that an interval is a point or degenerate interval if X= [x, x]. Also 
note that the left endpoint of an interior interval is always less than or equal to the 
right endpoint. The set of extended real numbers, i?* is the set of real numbers, R, 
extended with the two ideal points negative infinity and positive infinity: 

In the equations that appear in FIG. 5, the up arrows and down arrows 
indicate the direction of rounding in the next and subsequent operations. Directed 
rounding (up or down) is applied if the result of a floating-point operation is not 
machine-representable . 

The addition operation X+ 7 adds the left endpoint of Xto the left 
endpoint of Y and rounds down to the nearest floating-point number to produce a 
resulting left endpoint, and adds the right endpoint of Xto the right endpoint of Y 
and rounds up to the nearest floating-point number to produce a resulting right 
endpoint. 

Similarly, the subtraction operation Z- 7 subtracts the right endpoint of Y 
fi-om the left endpoint of X and rounds down to produce a resulting left endpoint, 
and subtracts the left endpoint of 7 from the right endpoint of X and rounds up to 
produce a resulting right endpoint. 

The multiplication operation selects the minimum value of four different 
terms (rounded down) to produce the resulting left endpoint. These terms are: the 
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left endpoint of X multiplied by the left endpoint of 7; the left endpoint ofX 
multiplied by the right endpoint of 7; the right endpoint of Z multiplied by the left 
endpoint of 7; and the right endpoint of X multiplied by the right endpoint of 7 
This multiplication operation additionally selects the maximum of the same four 
terms (rounded up) to produce the resulting right endpoint. 

Similarly, the division operation selects the minimum of four different 
terms (rounded down) to produce the resulting left endpoint. These terms are: the 
left endpoint of X divided by the left endpoint of 7; the left endpoint of X divided 
by the right endpoint of 7; the right endpoint of X divided by the left endpoint of 
7; and the right endpoint of X divided by the right endpoint of 7 This division 
operation additionally selects the maximum of the same four terms (rounded up) 
to produce the resulting right endpoint. For the special case where the interval 7 
includes zero, X/7is an exterior interval that is nevertheless contained in the 
interval 7?*. 

Note that the result of any of these interval operations is the empty interval 
if either of the intervals, X or 7, are the empty interval. Also note, that in one 
embodiment of the present invention, extended interval operations never cause 
undefined outcomes, which are referred to as "exceptions" in the IEEE 754 
standard. 

Bounding the Solution Set of a System of Linear Equations 

FIG. 6 illustrates the process of bounding the solution set of a system of 
linear equations in accordance with an embodiment of the present invention. The 
system starts by receiving a representation of the system of linear equations 
Ax = b, wherein A is an interval matrix and b is an interval vector (step 606). 

The system then preconditions Ax = b to produce Mqx = r, where 
Mo = BA and r = Bb (step 608). The preconditioning matrix B can be formed by 
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first computing an approximate center Ac of the matrix A, and then forming 
the approximate inverse of Ac, B = (Ac)"V Note that we can write 
A = Ac + Q[-l,l] where Q is a real matrix. Therefore, the preconditioned matrix 
is M = B A = I + BQ[- 1 , 1 ] . That is, the center of M is the identity matrix . If Ac 
5 and B were computed exactly, the center of Mo would be the identity matrix I. 
(Note that the system widens Mo at step 612 below so that the center of the result, 
M, equals I.) 

If we denote M = [mf{M),sup{M)] and r = [mf(r), sup{r)], then 
for ij = 1, infiMy) = -sup(My) (i j), and inf{Mu) + sup{Mu) = 2. 
10 In the preceding discussion, we have ignored the fact that B does not exist 

if Ac is singular. Suppose we try to invert Ac using Gaussian elimination. If Ac 
is singular, this fails because a pivot element is zero. At this point, the system 
terminates. Otherwise, if Ac is not singular and B can be computed, the system 
continues, 

1 5 Next, the system assures that supir^) ^ 0 by changing the sign of and Xt if 

necessary (step 61 1). Suppose we multiply the i-th equation of the system by -1 
and simultaneously change the sign of x,. As noted above inf(Mr) = -sup(My) for 
i ?t j. Hence, the off-diagonal elements are unchanged. Moreover, the diagonal 
elements change sign twice so they have no net change. Thus, the coefficient 
20 matrix is unchanged while x, and r, change sign. 

We can assure that supin) ^ 0 by changing the sign of and if 
necessary. Assume this is the case. If 0 g r„ we can change the sign of and x, if 
necessary and obtain -infin) < sup{r;). Therefore, we can always assure that 
0 < \inf(r,)\ < supin). Hereafter, we assume that the above relationship is satisfied 
25 for all /=!,..., n. This simplifies the procedure for finding the hull of Mqx = r. 

Next, the system widens Mo so that the center of the result, M, equals I 
(step 612). At this point, the system determines if M is regular (step 613). 
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If so, the system forms the hull h from M and r. In doing so, the system 
computes P= (M^)'^ as the inverse of the left endpoint of M (step 614). The 
system then forms = \l{2Pn - 1) for / = 1, . . (step 616), and also forms 
Zi = (mf{n)+sup{r;))Pn - e;^?sup{r) for / = 1 , . . . , wherein Cj^ is a unit vector in 
5 which the /-th element is 1 and other elements are 0 (step 618). Next, the system 
forms h by setting inf{ht) = c,z, if 2^ > 0 for = 1 , . . . , w, and by setting irif(k) = 
for / = 1, . . /? if < 0 (step 620). The system also sets supQn) = ?sup(r) 
(step 622). 

If M was not regular at step 613, the system uses the Gauss-Seidel process 
1 0 to compute the hull h (step 6 1 5) before terminating. 

The above-described procedure for finding the hull is valid only if M is 
regular. The following theorem enables us to verify regularity as a by-product of 
the computation of the hull 

15 Theorem 1 : Assume inf(M) is nonsingular so that P == (inf(M))"^ exists. Also 

assume that mXMO > 0 for all ? = 1, ... , n. Then M is regular if and only if P > I. 

If, using Theorem 1 , we find that M is regular, we can compute the hull h 
using the above-described procedure. Note, however, that the hull of the 
20 preconditioned system Mx = r is generally larger than that of the original system 
Ax = b. 

The foregoing descriptions of embodiments of the present invention have 
been presented only for purposes of illustration and description. They are not 
intended to be exhaustive or to limit the present invention to the forms disclosed. 
25 Accordingly, many modifications and variations will be apparent to practitioners 
skilled in the art. Additionally, the above disclosure is not intended to limit the 
present invention. The scope of the present invention is defined by the appended 
claims. 
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